library(Seurat)
library(tidyseurat)
library(tidyverse)
library(ggpubr)
library(patchwork)

sobj_sg <- read_rds('CRC-I/data/seekgene/crc-starsolo-annotated.rds')
sobj_10x <- read_rds('CRC-I/data/zy2020_tumor10x.rds')
sobj_sm2020 <- read_rds('CRC-I/data/smart_tpm.rds')
sobj_sm2018 <- read_rds('CRC-I/data/zhang2018_Count.rds')

# examine FCGR2B expression -------
sobj_sg |>
  VlnPlot('FCGR2B', group.by = 'zhang2020_fine', pt.size = 0, split.by = 'zhang2020_main') +
  coord_flip()

sobj_10x |>
  VlnPlot('FCGR2B', group.by = 'Sub_Cluster', pt.size = 0, split.by = 'Global_Cluster') +
  coord_flip()

sobj_sm2020 |>
  filter(!str_detect(Sub_Cluster, 'hE|hC|hF')) |>
  VlnPlot('FCGR2B', group.by = 'Sub_Cluster', pt.size = 0, split.by = 'Global_Cluster') +
  coord_flip()

sobj_sm2020 |> DotPlot(features = 'FCGR2B', group.by = 'Sub_Cluster')

sobj_sm2018 |>
  VlnPlot('FCGR2B', group.by = 'majorCluster', pt.size = 0) +
  coord_flip()

## FCGR2B in smart-seq2 -----------

### 2018 ------
fcgr2b_expr <- sobj_sm2018 |>
  get_abundance_sc_long('FCGR2B')

fcgr2b_expr <- sobj_sm2018 |>
  as_tibble() |>
  select(.cell, majorCluster) |>
  left_join(fcgr2b_expr)

fcgr2b_expr <- fcgr2b_expr |>
  mutate(label.main = str_extract(majorCluster, 'CD4|CD8')) |>
  filter(!is.na(label.main))

fcgr2b_expr |>
  ggplot(aes(label.main, .abundance_RNA, color = label.main)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type') +
  NoLegend()

g1 <- last_plot()

fcgr2b_expr |>
  group_by(label.main) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  count(label.main, wt = wt) |>
  ggplot(aes(label.main, n, fill = label.main)) +
  geom_col() +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells', color = 'Cell type') +
  NoLegend()

g2 <- last_plot()

g1 + g2 + plot_annotation(title = 'CRC TME T cell expression of FCGR2B (Zhang, 2018, Nature)')

#### sub CD8 ----------
fcgr2b_expr.cd8 <- fcgr2b_expr |>
  filter(label.main == 'CD8') |>
  mutate(label.fine = str_remove(majorCluster, 'h.+_|_C\\d{2}')) |>
  left_join(annotated_t_type)

fcgr2b_expr.cd8 |>
  filter(!str_detect(newtype, 'MAIT|Tn|GPR')) |>
  ggplot(aes(newtype, .abundance_RNA, color = newtype)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type') +
  NoLegend() +
  coord_flip() +
  scale_color_brewer(type = 'qual')

g1 <- last_plot()

fcgr2b_expr.cd8 |>
  group_by(newtype) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  count(newtype, wt = wt, .drop = FALSE) |>
  ggplot(aes(newtype, n, fill = newtype)) +
  geom_col() +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells in CD8 T cells', color = 'Cell type') +
  NoLegend() +
  coord_flip() +
  scale_fill_brewer(type = 'qual')

g2 <- last_plot()

g1 + g2 + plot_annotation(title = 'CRC TME CD8 T cell expression of FCGR2B (Zhang, 2018, Nature)')

### 2020 --------
fcgr2b_expr <- sobj_sm2020 |>
  get_abundance_sc_long('FCGR2B')

fcgr2b_expr <- sobj_sm2020 |> as_tibble() |>
  select(.cell, Sub_Cluster, Global_Cluster) |>
  left_join(fcgr2b_expr) |>
  filter(str_detect(Global_Cluster, 'CD4|CD8|Mye|B cell|ILC'))

fcgr2b_expr |>
  ggplot(aes(Global_Cluster, .abundance_RNA, color = Global_Cluster)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type') +
  NoLegend()

g1 <- last_plot()

fcgr2b_expr |>
  group_by(Global_Cluster) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  count(Global_Cluster, wt = wt) |>
  ggplot(aes(Global_Cluster, n, fill = Global_Cluster)) +
  geom_col() +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells', color = 'Cell type') +
  NoLegend()

g2 <- last_plot()

g1 + g2 + plot_annotation(title = 'CRC TME T cell expression of FCGR2B (Zhang, 2020, Cell)')

#### sub CD8 ----------
fcgr2b_expr.cd8 <- fcgr2b_expr |>
  filter(str_detect(Global_Cluster, 'CD8')) |>
  mutate(label.fine = str_remove(Sub_Cluster, 'h.+_|_C\\d{2}')) |>
  left_join(annotated_t_type)

fcgr2b_expr.cd8 |>
  filter(!str_detect(newtype, 'MAIT|Tn|GPR')) |>
  ggplot(aes(newtype, .abundance_RNA, color = newtype)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type') +
  NoLegend() +
  coord_flip() +
  scale_color_brewer(type = 'qual')

g1 <- last_plot()

fcgr2b_expr.cd8 |>
  filter(!str_detect(newtype, 'MAIT|Tn|GPR')) |>
  group_by(newtype) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  count(newtype, wt = wt, .drop = FALSE) |>
  ggplot(aes(newtype, n, fill = newtype)) +
  geom_col() +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells in CD8 T cells', color = 'Cell type') +
  NoLegend() +
  coord_flip() +
  scale_fill_brewer(type = 'qual', palette = 2)

g2 <- last_plot()

g1 + g2 + plot_annotation(title = 'CRC TME CD8 T cell expression of FCGR2B (Zhang, 2020, Cell)')

## in 10x -------
fcgr2b_expr <- sobj_10x |>
  get_abundance_sc_long('FCGR2B')

fcgr2b_expr <- sobj_10x |> as_tibble() |>
  select(.cell, Sub_Cluster, Global_Cluster) |>
  left_join(fcgr2b_expr) |>
  filter(str_detect(Global_Cluster, 'CD4|CD8|Mye|B cell|ILC'))

fcgr2b_expr |>
  ggplot(aes(Global_Cluster, .abundance_RNA, color = Global_Cluster)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type') +
  NoLegend()

g1 <- last_plot()

fcgr2b_expr |>
  group_by(Global_Cluster) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  count(Global_Cluster, wt = wt) |>
  ggplot(aes(Global_Cluster, n, fill = Global_Cluster)) +
  geom_col() +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells', color = 'Cell type') +
  NoLegend()

g2 <- last_plot()

g1 + g2 + plot_annotation(title = 'CRC TME T cell expression of FCGR2B (Zhang, 2021, Cell)')

#### sub CD8 ----------
fcgr2b_expr.cd8 <- fcgr2b_expr |>
  filter(str_detect(Global_Cluster, 'CD8')) |>
  mutate(label.fine = str_remove(Sub_Cluster, 'h.+_|_C\\d{2}')) |>
  left_join(annotated_t_type)

fcgr2b_expr.cd8 |>
  filter(!str_detect(newtype, 'MAIT|Tn|GPR')) |>
  ggplot(aes(newtype, .abundance_RNA, color = newtype)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type') +
  NoLegend() +
  coord_flip() +
  scale_color_brewer(type = 'qual')

g1 <- last_plot()

fcgr2b_expr.cd8 |>
  filter(!str_detect(newtype, 'MAIT|Tn|GPR')) |>
  group_by(newtype) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  count(newtype, wt = wt, .drop = FALSE) |>
  ggplot(aes(newtype, n, fill = newtype)) +
  geom_col() +
  theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells in CD8 T cells', color = 'Cell type') +
  NoLegend() +
  coord_flip() +
  scale_fill_brewer(type = 'qual')

g2 <- last_plot()

g1 + g2 + plot_annotation(title = 'CRC TME CD8 T cell expression of FCGR2B (Zhang, 2021, Cell)')


